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ABSTRACT 



A computer program was developed to solve the quasi-one- 
dimensional, unsteady flow problem associated with starting 
phenomena in various nozzle-short diffuser combinations of 
interest to gas dynamic lasers. The program was based on the 
techniques developed by George Rudinger in his book, Wave Dia - 
grams for Nonsteady Flow in Ducts . The techniques were modi- 
fied to facilitate a fixed grid for computational ease. 

The program was used to analyze the starting transients 
of a Mach 3 nozzle with a semi-wedge diffuser. Results indi- 
cated that for a given inlet /exhaust pressure ratio, diffuser 
starts could be obtained at higher ramp angles . for thin dif- 
fusers than for thick diffusers. 
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I . INTRODUCTION 



Gas Dynamic Lasers (GDLs) typically operate at cavity 
Mach numbers of 4.0 to 5.5 in order to obtain the kinetics 
desired for population inversion and signal gain. GDL opera- 
tion often consist of a series of short bursts which result 
in periods of supersonic flow followed by periods of no flow. 
Although efficient pressure recovery is not of prime impor- 
tance, the high Mach number, low pressure flow must be dif- 
fused in order to permit exhausting to the atmosphere. While 
rapid and efficient diffuser start-up is required in all 
laser applications due to its effect on lasing performance, 
it becomes especially critical in the short burst mode of 
operation. 

Some chemical lasers also require diffusing of low pres- 
sure high flow rate gases. Since these gases are collected 
and condensed or absorbed, diffuser efficiencies become more 
important in this type of application. 

Aircraft installation of lasers compounds the above prob- 
lems with size, weight and pressure ratio limitations. Dif- 
fusers used on GDLs currently constitute a large percentage 
of the laser size and hence become prohibitive in airborne 
installations. If an aircraft axial flow compressor is used 
for gas supply, pressure ratios are limited to approximately 
30 to 1. 

The diffuser starting problem is similar to that obtained 
in starting supersonic wind tunnels. Extensive work was done 
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in this area following World War II and into the early 1950 's. 
Some of the means which have been utilized in starting these 
tunnels are variable area diffusers, overpressure, perforated 
walls, auxiliary air injection, and boundary layer control 
[Ref. 6]. 

In conjunction with work being done for NAVAIR at China 
Lake on the airborne gas dynamic laser, research is being 
undertaken at the Naval Postgraduate School on minimizing 
diffuser size requirements and starting time. The research 
has proceeded along two directions; experimental and analyti- 
cal. The experimental investigation of diffuser start-up has 
been conducted by LCDR J. Zerr [Ref. 8], utilizing a blow- 
down wind tunnel to establish Mach 3 and Mach 4 flow. As a 
first step in the analytical program it was desired to have 
a computer model by which various flow conditions could be 
analyzed. 

A complete analysis of the starting process would require 
a four-dimensional approach; 3 spatial coordinates plus the 
time coordinate. In such an approach, discontinuities and 
non-linearities caused by the presence of shock waves create 
problems of enormous complexity. 

One can obtain insight into the starting phenomenon through 
a quasi-one-dimensional analysis which utilizes one space and 
one time coordinate. Transverse area-change effects may be 
included if the radius of curvature of the change is large in 
comparison with the axial duct length. 
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Rudinger [Ref. 5] treats the one-dimensional unsteady, 
non-linear wave problem for flow in ducts. The procedures 
he develops allow for the inclusion of heat addition, vari- 
able entropy, variable area and several boundary conditions 
and discontinuities. Since his method is designed for hand 
graphical analysis, however, it does not lend itself directly 
to computer numerical methods. 

Spalding [Ref. 7] has modified the techniques of Rudinger 
in order to permit efficient digital computation procedures. 
His modifications consist primarily of the implementation of 
a fixed spatial grid combined with a backward interpolation 
of wave characteristics. Additionally, the Riemann variables 
for the characteristics have been re-defined for ease of com- 
putation . 

Solution of the characteristic-wave problem is also dis- 
cussed extensively in Refs. 1, 3, and 6. Several alternatives 
are offered for handling the iteration process which inevit- 
ably results from the inclusion of the area effect term. The 
complexity of these schemes, in any case, must be weighed 
against the increase in accuracy and the obvious limitations 
of the one-dimensional analysis. 



13 



II. APPROACH 



The computer program developed in this thesis is based on 
the quasi-one-dimensional unsteady wave analysis described in 
Ref. 5 as modified by Ref. 7. These techniques have been 
utilized in the development of a previous computer program, 
PREX, for use with cyclic, constant area pressure exchangers. 
Since many of the subroutines used in PREX were general in 
nature, they were adapted as a starting point for this thesis 
and extensions were developed to permit the calculation of 
supersonic flow through ducts of varying area. 

The program essentially solves a shock-tube type problem 
in a duct with one closed end and one open end, (i.e., a 
Ludweig tube). The nozzle and diffuser are located between 
the diaphragm and the open end. When the diaphragm is broken, 
pressure waves travel through the nozzle/diffuser section 
while rarefaction waves proceed toward the closed end. Since 
the diaphragm is located closer to the open end, the diffuser 
starts prior to the time the rarefaction waves reach the 
closed end. In this manner the diffuser start-up process is 
isolated from the reflected rarefaction waves returning from 
the closed end in order to more closely approximate a reser- 
voir boundary condition. The computational model is depicted 
in Figure 1. 
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Figure 1. Duct Geometry. 



III. DERIVATION OF EQUATIONS 



In the derivations of the characteristic equations, the 
following implicit assumptions have been made; 

1. The gas flow obeys the ideal gas relations. 

2. The gas is calorically perfect (constant specific 
heats) . 

3. The flow is considered quasi-one-dimensional (area 
effects are included). 

4. The duct is rigid in time with no mass flux through 
the walls (e.g., no blowing or suction). 



A. CONTINUITY 



Figure 2 depicts the mass entering and leaving a control 
volvime during time dt. Conservation of mass requires: 

Increase of mass within Net mass flux through 

the control volume the control surface 



•^(pAdx)dt = pAudt-[pAu + -^(pAu)dx]dt 

o L dX 

|j(pA) + I^(PAU) - 0 



Since the duct is rigid, A is not a function of time but 
is a function only of x. Therefore expansion of the above 
equations- results in: 



i£ + PUdA 
9t Adx 





0 
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pAudt 

+ -|^(pAu)dxdt 



Figure 2, Mass Flow Through Control Volvime During Time dt. 




Figure 3. Forces Exerted on Control Volume 
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B. MOMENTUM 



Figure 3 depicts the forces on a control surface within 
a duct of varying area. From Newton's second law: 

body acceleration = acting on the body 

mass of the body 

Let f = summation of frictional dissipative forces/unit mass; 
then 

Expanding the above equation, one gets 

li ^ “IS “ -sidsip* ^ * p#>’‘ 

iu + uHi = _ilE+f 
at ax p ax 



C . ENERGY 

The preceding control volume approach may also be used 
for a similar derivation of the applicable energy equation. 
Several reliable references have used this approach and the 
results will not be duplicated here; exceptionally thorough 
discussions are contained in Ref. 3 and 8. For an adiabatic, 
quasi-one-dimensional flow of a perfect gas the energy equa- 
tion reduces to: 

i_(e + ^) + u-^(e + Hi) + i llPJil + EH_dA ^ q 
9t^ 2 ' ^3x^ 2 ' p 9x Ap dx 

where e is the gas internal energy/unit mass. 



D. CHARACTERISTIC EQUATIONS 

The equations of continuity, momentum and energy may be 
combined to obtain the characteristic (Riemann) wave equations. 
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Multiplication of the momentum equation by u yields 



8u , , u 9p o 

u-:^ + u^tt- + — = fu 

8t 8x p 9x 

Subtracting the above from the energy equation, and assuming 

that f=0, gives 

9e ^ ^ P ^ pu dA _ , 

9t 9x p 9x pA dx 

The continuity equation 

i£ + dA ^ 9u ^ ^ 9£ ^ Q 



9t 



A dx 



9x 



9x 



may be re-written in the form 

p9£^£U9p^£^^£U^^ 
p 9t p^ 9x p 9x pA dx 

Subtracting this from the energy equation yields 

+ Je _Pl£_£u9£=n 
9t 9x ^ 9t ^ 8x 

For an ideal gas with constant specific heats 

,/S. 2 da dp 

^<R> ' 0 ■ 

For quasi-one-dimensional adiabatic flows the entropy of each 
particle may be expressed as 



Combination with the continuity- equation produces 

If * “H * • 

The momentum equation 

9u + u^ + iiE=0 

9t 9x p 9x 

combined with the perfect gas relation 



( 1 ) 



dp _ 2y da _ ^ 
p Y~1 a R 
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yields 



+ u^) + a— - 9(s/R) = 0 

2 ' 8t 9x"^ 9x 2y 9x 

or 

Y-1 ^ 9u . SUv . „9a _ a^ 9s 

2 <H ^-35 - 2S" -SS • 

Jb*^ 

Addition and subtraction of Equations 1 and 2 gives 



( 2 ) 



a‘ 



|j(a + u) + (ufa)|j(a + ^ u) = 

9 / Y “1 \ j_ / \ 9 /" ^“1 \ _ 3 -‘ 

3t^^ 2 ^^ 9 x^^ 2 2 c_ 9 x 



9s _ ( Y~l )ua dA 



dx 



3s _ ( Y~l )ua- dA 
2 A dx 



Note that the above two equations can be interpreted as the 

rate of change with time of the quantities (a + and 

Y-1 1 1 

(a - -hr- u) along lines with slopes of -r- and respectively 

^ VI * VI 

in the x, t plane. Defining new variables [Ref. 7] 



U = ^u 



P = (p) 



a = exp [(s-s^)/ 2 Cp] 



where 



s* = So + 2c In ao - R In P( 



♦ -p 

which implies 
Po = a 

the derived equations may be re-written 

Ij(Pc-hO) + (u+a)|3;(Pa-^0) = |h |f ‘ f 

9 / r > TT \ ^ \9 /- r > TTN a^ 3s ^ Y-l \Ua dA 

-^(Pa-U) + (u-a)^(Pa-U) = - 7 ^ - i-^)~ 



3x 



2Cp 3x 



A dx 
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Referring to Figure 4 and applying the above equations it 
can be seen 



P a +U 
P P P 


' '’qV"q * 


•^PQ 


a2 

2“p 


as 

ax 


^PQ 


2 ’k 


dx 


P a -U 
s s s 


= P a -U - 
s s s 


'ps 


2"p 


as 

ax 


•^ps 


2 ’k 


dA 

dx 



which are the Riemann equations along the characteristics ex- 
pressed in terms of the newly defined variables P , a , and U . 
Reference 2 treats the constant area case of duct flow for 
which the above equations become 



P 0 +U 




/ 




as 


P P P 


Q Q Q 


pQ 


2c 

P 


ax 


P 0 -U 


= P 0 -U - 


/ 




Is 


s s s 


s s s 


ps 


2c 


ax 



Through an analysis of wave interaction including a con- 
tact discontinuity, a justification is presented for writing 
the above equations as ' ^ 



P Q P Q Q Q 

P a -U = P a -U 
p s p s s s 



and Op = that value of a determined from the particle path of 
slope 1/u. 

Let HPS and HPQ represent the average values of 

^ along the characteristics PS and PQ respectively. 
Then HPS dt and HPQ dt represent changes along the character- 
istics over time step dt due to the changing area. The term 
^ may be easily obtained at points S and Q and an 
iterative procedure allows it to be evaluated at point P. 
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Figure 4. Computational Procedure. 
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The equations may now be re-written as 



P 0 ^+U = +U^-HPQ dt 

P Q P Q Q Q 



P 0 -U = P 0 -U -HPS dt 
p s p s s s 



which may be solved for P and U as 

P P 

P^ 0 ^+P 0 +U^-U -HPQdt+HPSdt 
p ^ Q Q s s Q s 

p ■ (V°s> 



( 3 ) 



.. PQ-P3+UQ/PQ+U3/«3-HPQdt+HPSdt 
P (l/Oq + I/O 3 ) 



(4) 



0 = 
P 



R 



(5) 



where R is the x-intercept point of a characteristic having 
a slope 1 /u. 

Equations (3), (4), and (5) are the equations which are 
actually programmed. They are utilized in subroutine STEP 
and are modified by boundary conditions for use in determin- 
ing duct end values in subroutine PORTS. A discussion ' of 
these modifications as well as how the equations are utilized 
is contained in a description of the program. 
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IV. DESCRIPTION OF PROGRAM SS DIFFUSER 



A. PROGRAM EVOLUTION 

The computer program developed is based on a program 
called PREX by Lebius Matthews [Ref. 7]. PREX (mnemonic for 
pressure exchanger) solves the unsteady flow problem result- 
ing from opening and closing several valves which separate 
the duct from ports maintained at various pressure levels. 

The valves utilized are modeled as orifices which change 
in size linearly with time during the opening and closing 
process. Extensive work was done in attempting to adapt this 
concept for application in program SS DIFFUSER (mnemonic for 
super-sonic diffuser). In this application extremely rapid 
valve operation was desired in order to simulate the high 
pressure rises obtained in laser starting. Additionally, 
the valves had to be re-modeled to handle subsonic and super- 
sonic flow. 

In all our attempts to use the valves in PREX, erroneous 
stagnation pressure increases were observed near the area of 
the valve. Once established, these new stagnation pressures 
then propagated into the nozzle/diffuser section. Addition- 
ally the magnitude of the established flow Mach number in the 
duct seemed to be a function of the valve opening time. The 
valves were finally abandoned in favor of a more simplified 
boundary condition which consists of a bursting diaphragm 
such as witnessed in shock tube applications. Subroutine 
PORTS was subsequently completely re-written to handle the 



24 







I 



boundary conditions as outlined by Rudinger [Ref. 5]. Minor 
modifications were made in subroutine PORTX and in the main 
program, combined with the deletion of two other subroutines 
which handled the cyclic aspect of the pressure exchanger 
problem. 

Most of the changes made to handle the variable area duct 
were made in subroutine STEP. This is also the subroutine 
which would have to be changed to include frictional or heat- 
ing effects. 

Subroutine STEP utilizes an entropy averaging process to 
avoid the discontinuity problems associated with shock waves. 
This is developed and discussed in depth in Ref. 6. The im- 
pact of this averaging on the results is that entropy values 
near the shock can not be I'elied upon to produce an accurate 
representation of the physical process. Also, if a station- 
ary shock appears somewhere in the duct, the entropy gradient 
eventually averages out entirely. 

B. EXPLANATION OF SUBROUTINES 

The following is an explanation of subroutines used in 
program SS DIFFUSER. Figure 5 contains a flow chart of the 
computational sequence. A list of symbols used in the com- 
puter program can be found in Appendix B. 

1. BLOCK DATA (data input) 

BLOCK DATA are divided into six categories of input 
data. Under FLUID PROPERTIES, the specific heat, Cp and c^, 
of the fluid are specified. These are treated as constants 
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Figure 5» Flow Diagram of Program SS DIFFUSER. 
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throughout the program. Input data under BOUNDARY CONDITIONS 
are the temperatures and pressures at the ends of the duct. 
IDATA allows one to specify static or stagnation pressure in- 
put values. In the sample program these have been non-dimen- 
sionalized with respect to atmospheric conditions. In section 
VALVE PARAMETERS, the input value of NOVALV allows one to 
specify a fully open or fully closed right duct end. INITIAL 
CONDITIONS data fixes the values of temperature and pressure 
within the duct prior to the start of calculation. 

DUCT PHYSICAL GEOMETRY information is used by subroutine 
ARAT to determine the area effect on the flow. A more de- 
tailed discussion of these values is included in the descrip- 
tion of subroutine ARAT. PLOTTING AND PRINTING PARAMETERS 
controls the printing and plotting frequency and spacing of 
spatial nodes in the output. It also provides for a maximum 
time of solution, TMAX. 

2 . MAIN (main program) 

The main program controls the order of calculations 
and specifically the calling of the designated subroutines. 

All the calculations are done on the subroutine level and 
passed to main through dummy variables and common statements. 

MAIN utilizes TMAX to determine the maximum solution 
time. / 

3 . INIT (initialization) 

This subroutine utilizes the specific heat inputs to 
determine the GAM- constants (see Program Terminology, Appen- 
dix B) used throughout the program. It also establishes the 
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spatial interval according to the number of calculation points 
desired. 



Slopes, areas, temperatures and pressures are initial- 
ized to the input values. From this information, initial 
values of the entropy variable, SIG, are determined. 

4. ARAT (area ratio) 

Subroutine ARAT utilizes the duct geometry to calcu- 
late the value of the area effect term in the Riemann Equa- 
tions. Its output, RAT, is a ratio of the slope of the duct 
wall to the duct area at each value of x. Areas and x-coor- 
dinates have been non-dimensionalized with respect to duct 
length in the sample program. 

To provide for a smooth change of area in the nozzle 
section, the area change has been modeled as a sinusoidal 
variation. The sinusoid amplitude, B, is equal to one-half 
the difference between the original duct area and the area at 
the nozzle throat. Thus, a smooth contour having zero slope 
at nozzle entrance, throat, and exit is provided. 

The diffuser section allows for the modeling of a dif- 
fuser having an initial ramp input and three subsequent changes 
in flow direction. Figure 1 depicts the above models. 

5 . AUX (auxiliary) 

Subroutine AUX calculates the values of output quan- 
tities other than at spatial calculation nodes. Although the 
values of A(I) and U(I) are used by subroutine STEP, other 
computations made here are not critical to calculations made 
by other routines. Their inclusion serves primarily as a 
check on input data and boundary conditions. 
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6 . OUTPUT (printed output) 



This routine gathers information from other subrou- 
tines and consolidates it for printing and plotting purposes. 
The plotting routine, UTPLOT, is a catalogued routine avail- 
able under the IBM 360 system which was utilized. Similar 
such routines may be easily substituted. 

The variables shown in the sample output have been 
non-dimensionalized with respect to atmospheric conditions. 

The entropy value (s) is actually s/R. 

7. STEP (time step) 

In a sense, STEP is the real controlling routine with- 
in the program. It controls the calculation at all interior 
points and uses routines SLOPE and LININT to do so. It es- 
sentially uses a backward linear interpolation scheme which 
may best be described by reference to Fig. 4. 

Consider the calculation at some interior point P. 

At some previous time, ti, the solution exists at points Q, 

R, and S as well as along the entire spatial grid. Since it 
is desired to maintain a constant grid interval, character- 
istics 1, m, and n are generated through P based on the pre- 
viously known values of U and A at points Q, R, and S. For 
the m characteristic shown, this would be the average of the 
values of U+A at Q and R. This same techniques is then ap- 
plied along the entire length of duct. 

The linear interpolation scheme requires that these 
characteristics intersect the x-axis at time ti somewhere be- 
tween the calculation node and the adjacent node. To 
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facilitate this, the time step is adjusted to correspond with 
the minimiun slope characteristic. Additionally, since this 
characteristic is only a first approximation, only a portion 
of this step is utilized as dictated by duct geometry and the 
expectation of area effect. It may be readily seen that the 
calculation time is quite sensitive to the number of grid 
points . 

Once the time interval has been established, the cor- 
responding x-intercepts are calculated by subroutine SLOPE. 
Values of the computing variables are then calculated at these 
intercepts by linear interpolation of the bracketing calcula- 
tion nodes Q, R, and S. Using these values, a rough guess 
for the calculation variables are obtained for point P based 
on the Riemann equations with no area change. 

The area effect term is integrated along the charac- 
teristic by assuming it to be a constant equal to the average 
of the end values. Subsequently new values are obtained at 
P through the introduction of this term. An iteration pro- 
cedure is introduced here in order to obtain a more accurate 
area effect term. 

Upon completion of the velocity and pressure calcula- 
tion at point P, the entropy term, SIG, is updated at time 
t + dt . 

8. SLOPE (characteristic slopes) 

As previously mentioned, subroutine SLOPE is respon- 
sible for computing the slopes of the characteristics. Char- 
acteristics which pass through endpoints of the duct are 
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calculated only for the one originating within the duct. 

These are then combined with boundary values in subroutine 
PORTS (Fig. 4). 

Two options may be used in calling SLOPE. The first 
option computes a slope on the basis of the average values 
of the bracketing spatial nodes. This is the option which is 
utilized in the sample program. The other option, which may 
be used with appropriate changes in subroutine STEP, allows 
the values in the slope equation to be externally dictated. 

In either of the above options, the subsequent x- 
intercept at time t-dt is computed based on the slope and 
time interval. 

9. LININT (linear interpolation) 

LININT performs a simple linear interpolation for 
variables based on the values of the x-calculation nodes 
bracketing the desired spatial coordinate. The interpolation 
always takes place between two adjacent nodes. Endpoint ex- 
clusion is accomplished in the same manner as described for 
subroutine SLOPE. 

10 . PORTX (port auxiliary) 

PORTX is the controlling routine for calculation of 
pressure and temperature variables external to the duct. If 
the duct ends are closed, the dummy variables passed to PORTS 
are appropriately set to zero. If the duct end is open, as 
in the case of the sample program, the variables are calcu- 
lated based on whether static or stagnation values have been 
specified. These values arc then passed to PORTS where they 
are meshed with the characteristic values. 
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11 . PORTS (port solutions) 



Subroutine PORTS controls the calculation of the ve- 
locity and pressure variables at the duct ends. If the ends 
are closed, PORTS provides the appropriate wave reflection 
from the Riemann equations. 

For outflow, PORTS differentiates between two cases, 
sub-critical (subsonic) and super-critical (supersonic). For 
the sub-critical case the existing outside pressure is uti- 
lized as discussed in Ref. 5. For supersonic flow, all char- 
acteristics originate from within the duct and hence the 
Riemann equations are directly applicable. PORTS also pro- 
vides an indicator to OUTPUT indicating whether the flow is 
supercritical or subcritical. 
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V. DISCUSSION OF RESULTS 



The computer program was initially used to solve a sample 
problem from Ref. 5. In Appendix A the computed solution is 
presented along with the graphical solution given by Ref. 5. 
The results can be seen to be in very good agreement with 
those of Rudinger. 

Additional flows were considered utilizing nozzle con- 
traction ratios of 4:1 and gas specific heat values c^ and c^ 
of 0.172 and 0.24 BTU/lbm-°R, respectively. Isentropic values 
for these conditions indicate a resulting Mach number of 2.94 
should be obtained. Mach numbers of approximately 3.04 were 
obtained indicating an error of 3.4%. This resulting error 
is quite low considering the inherent limitations of the lin- 
ear interpolation process utilized. 

Tables I-V depict typical results obtained for a driver 

pressure value below that required to obtain a fully started 

a„t ^ 

condition. Time was nondimensionalized ( “ ) utilizing a four 

Jjo 

foot duct length. The pressure utilized was p/po of 16.0 
(Po = ambient pressure). A semi-wedge diffuser (see Fig. 6) 
was utilized having a ramp slope of 0.05 (6=3®) and contrac- 
tion ratio ( ) of 0.94. Subsonic flow is es- 

^test section 

tablished ahead of the nozzle throat (x=.8125) with supersonic 
flow in the diverging section. A formation of a stationary 
shock wave (Fig. 6) can be observed in the diverging section 
(x=.8700) of the nozzle resulting in subsonic flow in the 
test section (x=.8750 to x=.8900). 
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Tables VI-IX depict values for the same duct geometry for 
which the driver pressure has been elevated to 20.0 to drive 
the shock wave out of the nozzle and downstream (Fig. 7). 

Pressures necessary for starting were considerably greater 
than witnessed in the laboratory. Possible explanations for 
this are that the program incorporates no friction, the burst- 
ing diaphragm causes a stagnation pressure loss, and only nor- 
mal shocks are modeled. An effective stagnation pressure 
ratio of approximately 18.0 can be observed in the nozzle. 

When combined with the losses across a Mach 3.0 shock 
(Pt 2 /Pti“* > a correlation can be seen with laboratory ob- 
served starting values of 6.0. 

These results appear to be in good agreement with the 
cited references. References 1 and 4 discuss similar condi- 
tions extensively. Although they treat the diverging section 
as a step increase in duct area, similar results are reported. 
For weak shock waves passing through a widening duct, a re- 
flected expansion wave is created along with a transmitted 
shock. Additionally, a contact discontinuity is created by 
the entropy increase due to the shock. The resulting flow is 
subsonic in all regions; the entire process is shown in Fig- 
ure 8A. For stronger shock waves region 5 disappears (Figure 
8B), and at yet stronger values a stationary shock occurs as 
shown in Figure 8C. Further increases in arriving shock 
strength cause the stationary wave to be swept downstream as 
depicted in Figure 8D. Although these distinct waves can not 
be easily picked out in the computer results due to complica- 
tions caused by a continuously varying area and the effects 
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of the nozzle contraction preceding in the area of interest, 
the general trends described above are in qualitative agree- 
ment. It appears as though the area contraction has a strength- 
ening effect on the shock and causes an increase in shock 
velocity as might be expected from continuity. This even- 
tually results in the formation of a reflected wave travel- 
ling upstream. Such an effect can be seen when a diffuser is 
added to a nozzle with the resulting larger pressure gradient 
across the stationary shock wave in the nozzle diverging sec- 
tion . 

Figure 9 shows the effect of diffuser contraction ratio 
and ramp angle on start-up. For a fixed nozzle entrance to 
diffuser exit pressure ratio of 20.0, startups were attempted 
utilizing 56 combinations of ramp angle and contraction ratio. 
The results indicate that thin diffusers permit the use of 
higher ramp angles than do thick diffusers. Zerr [Ref, 8] 
determined similar results experimentally and advocates the 
use of thin, high ramp angle diffusers for efficient starting. 
Values for high ramp angles and extremely thin diffusers are 
not shown since these values slip through the computational 
grid size utilized (x=.01). 
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Table I. Example 1: Initial Conditions 



PRESSURE DENSITY VELOCITY TEMP ENTROPY MACH 
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Table II. Example 1: Conditions at Time = .0599 



PRESSURE DENSITY VELOCITY TEMP ENTROPY MACH 
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Table III. Example 1: Conditions at Time = .1438 



PRESSURE DENSITY VELOCITY TEMP ENTROPY MACH 
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Table, IV. Example 1: Conditions at Time = .1808 



PRESSURE DENSITY VELOCITY TEMP ENTROPY MACH 
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Table V. Example 1: Conditions at Time = .2364 
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Figure 6. Example 1: Pressure Distribution at Time - .2364 
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Table VI. Example 2: Conditions at Time = .0505 



PRESSURE DENSITY VELOCITY TEMP ENTROPY MACH 
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Table VII. Example 2: Conditions at Time = .0724 



PRESSURE DENSITY VELOCITY TEMP ENTROPY MACH 
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Table VIII. Example 2; Conditions at Time = .1926 



PRESSURE DENSITY VELOCITY TEMP ENTROPY MACH 
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Table IX. Example 2: Conditions at Time = .2243 
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Figure 7. Example 2: Pressure Distribution at Time = .2243 
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Figure 8. Waves Resulting From Shock Wave Passing Through 

a Widening of the Duct as a Function of Increasing 
Pressure Ratio. 
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Figure 9. Effects of Diffuser Ramp Angle and Contraction Ratio on Starting, 



VI. CONCLUSIONS AND RECOMMENDATIONS 



The developed computer model appears to give satisfactory 
results in one-dimensional wave analysis of the diffuser 
starting phenomena. Only a few of the possible parameters 
have been varied thus far but the program may be used for 
more extensive studies in a number of areas. 

The high contraction ratio nozzle did not appear to vio- 
late the implicit quasi-one-dimensional assumption. Reference 
1 gives the criterion that the axial component of the flow 
velocity and its derivatives should be larger than the cor- 
responding transverse component by at least one order of mag- 
nitude, to satisfy this assumption. While this criterion may 
be overly restrictive in some applications, it must be abided 
by in high contraction ratio problems. 

The program appears to have no difficulties in the tran- 
sonic range where one of the characteristics develops a ver- 
tical slope. A problem was anticipated in this regime, and 
slope and interpolation equations were written so as to avoid 
indeterminate forms. 

One of the most restrictive features of the computer pro- 
gram is the time step dependence on grid size. The interpo- 
lation process requires that the characteristics at time step 
t + dt originate within one spatial grid space at time t. 

Thus in decreasing the grid size for more accuracy one must 
pay a high price in calculating time. 
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As previously mentioned, extensions may be made in sub- 
routine STEP to include frictional and heat transfer effects. 
Additional changes may be possible to include other boundary 
layer phenomena. 

Studies so far have been based on diffusers with a lead- 
ing ramp angle followed by a constant area section to the duct 
end. A similar study may be done on a full wedge diffuser 
(converging/diverging) in order to gain a more complete pic- 
ture of diffuser effects. 

Reference 2 documents the results of several applications 
of the basic techniques utilized by the program. These in- 
clude a shock tube, supercharger and pressure exchanger. 
Program SS DIFFUSER may be used for constant area shock tube 
problems by setting the diffuser wedge angles and nozzle 
sinusoid equal to zero and utilizing the closed end condition 
for the right end. 
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APPENDIX A 



RUDINGER DUCT PROBLEM 

The following problem was taken from Ref. 4. The compu- 
ted output is presented for comparison with the graphical 
solution presented by Ref. 4. 

PROBLEM : 

A duct of constant cross section is filled with air 
(Y=1.4) that is isentropically compressed from atmospheric 
pressure to Pi =2. 83. The duct is suddenly opened to the at- 
mosphere. How does the pressure vary at the closed end of 
the duct? 

SOLUTION; 

The pressure history at the duct end is depicted in Fig- 
ure Al. Close agreement is seen throughout. Initially an 
expansion wave travels toward the closed end and is reflected 
back to the open end. The open end condition causes it to 

reflect back as a shock wave which reaches the closed end 
at 

shortly after ■ = 3.0. 

There is a slight variation of results in the vicinity of 
the shock wave where the computer solution depicts a finite 
slope instead of the predicted vertical discontinuity. This 
is a result of the averaging process across the shock previ- 
ously mentioned. 
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Figure Al. Rudinger Duct Problem Results 



APPENDIX B 





PROGRAM TERMINOLOGY 


Fortran Symbol 


Meaning 


A(I) 


a Sonic velocity 


AA 


Duct area 


AK 


Fraction of maximum allowable time step 


AM(I) 


U Velocity variable 


AMAC(I) 


Local Mach number=U/A 


AMI 


Initial 'AM' in duct 


AMl(I) 


'AM' interpolated along 'N' character- 
istic 


AM2(I) 


'AM' interpolated along 'M' character- 
istic 


AMT 


'AM' prior to inclusion of area effect 
term 


AN(I) 


Pressure variable 


^ ANl 


Initial value of 'AN' in duct 
(low pressure side) 


ANII 


Initial value of 'AN' in duct 
(high pressure side) 


ANl(I) 


'AN' interpolated along 'N' character- 
istic 


AN2(I) 


'AN' interpolated along 'M' character- 
istic 


AREA 


Duct area (in ARAT) 


ANT 


'AN' prior to inclusion of area effect 
term 


AUX 


Subroutine name 


B 


Amplitude of nozzle sinusoidal area 
variation 



N 
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BE(I) 




Characteristic slope 


BEM(I) 


B 

m 


Slope of 'M' characteristic 


BEN(I) 


B 

n 


Slope of 'N' characteristic 


BES(I) 


B 

s 


Slope of 'S' characteristic 


BMAX 




Maximum of BEM, BEN, and BES 


C 




Label for common block 


Cl 




Label for common block 


C2 




Label for common block 


C3 




Label for common block 


CELL 




Label for common block 


COMP 




Label for common block 


CP 


c 

P 


Specific heat at constant pressure 


CV 


c 

V 


Specific heat at constant volume 


CYC 




Label for common block 


DADX 




Da/dx change in duct area with respect 
to X 


DG 




Increment in 'G' 


DL 




Duct length 


DT 




Time step size 


DX 




Spatial step size 


EPS 




Iteration control term 


FLOW 




Label for common block 


FRA 




Fraction in linear interpolation 


G(l>IPORT) 




Integrated mass flow into port numbered 
' NPORT ' 


GAM 




Ratio of specific heats=CP/CV 


GAMC 




Critical pressure ratio for 
INFLOW=SQRT( (GAS+1 )/2) 
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GAMMA 


Label for common block 


GAMl 


= (GAM-l)/2 


GAM2 


= -2/GAMl 


GAMS 


= 2(GAM-1)/(GAM+1) 


GAM4 


= (2/(GAM+1))**((GAM+1)/(2*(GAM-1)) ) 


GAMS 


= GAM/ GAMl 


GAM6 


= 1/GAM5 


GAM7 


= (GAM+l)/(2+GAM) 


GAMS 


= (GAM+1)/ (GAM-1) 


HPS 


Area effect term along 'M' character- 
istic 


HPQ 


Area effect term along 'N' character- 
istic 


IDATA(I) 


INDEX^ = 1 If port static pressure supplied 

= 2 If port stagnation pressure supplied 


IHAND 


Index identifying port 1 for left hand 
port, and 2 for right 


IND 


Indicator to route computation 1 for NC 

flow closed port 8 for super-criti- 

cal outflow through fully open port 9 
for sub-critical outflow through fully 
open port 


INDL 


'IND* at left end 


INDR 


' IND' at right end 


INIT 


Subroutine name 


I OPEN 


Index, =2 if end closed, =1 if end open 


IP 


Alphameric symbol identifying end 


IR 


Alphameric symbol for right end 


J 


Counter 


K 


Counter in DO loop 


K1 


KDPM + 1 



KDPM 


Location of diaphragm in duct 


L 


Causes print of variables every L-TH 
node 


LININT 


Subroutine name 


MAIN 


The main program 


N 


Numbers of grid intervals in spatial 
direction 


N1 


Numbers of grid points in spatial 
direction = N+1 


NOPENL 


Number of port open on left hand side 


NOPENR 


Number of port open on right hand side 


NOVALV 


NOVALV = 1 Implies right end open, 
closed for all other values 


NPLOT 


Plot control 


NPRINT 


Print control 


NSTEP 


Number of time steps 


OUTPUT 


Subroutine name 


P(I) 


p Pressure 


PHI 


Duct end opening (if open, if closed) 


PHI 


Sinusoidal angle in ARAT 


PHIL 


'PHI' at left port 


PHIR 


'PHI' at right port 


PI 


Initial pressure in duct (low pressure 
side ) 


PII 


Initial pressure in duct (high pressure 
side) 


PIE 


3.1415926536 


PO 


Pg^ Stagnation pressure 


POOL 


Stagnation pressure in duct, left 
boundary 



POCR 


Stagnation pressure at duct, right 
boundary 


POL 


Stagnation pressure at left end (exter- 
nal) 


POP(I) 


Stagnation pressure at duct end (exter- 
nal) 


POR 


Stagnation pressure at right end 
(external) 


PORTS 


Subroutine name 


PORTX 


Subroutine name 


POS 


Non-dimensional spatial coordinate 


PR 


Ratio of stagnation to static pressure 


PRL 


'PR' at left end (external) 


PRR 


'PR' at right end (external) 


PROP 


Label for common block 


PX 


p^^ Static pressure (external) 


PXL 


Static pressure at left end (external) 


PXP(I) 


Static pressure in port 'I' 


PXR 


Static pressure at right end (external) 


PXS 


Dummy argument for 'PXL' or 'PXR' 


RAN(I) 


Limits of output plot 


RANGE 


Label for common block 


RATIO 


Label for common block 


RAT(I) 


Dummy variable in ARAT for ^ 


RATX(I) 


RAT evaluated at C(I) 


RATXl(I) 


RAT interpolated along 'N' characteristic 


RATX2 ( I ) 


RAT interpolated along 'M' characteristic 


RO(I) 


p Specific density 


S(I) 


s Entropy 
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S(I) 




Dummy variable in linear interpolation 


S1(I) 




Dummy variable in linear interpolation 


SIG(I) 




Dependent variable, function of entropy 


SIGl(I) 




SIG interpolated along 'N' characteristic 


SIG2(I) 




SIG interpolated along 'M' characteristic 


SIGI 




Initial value of SIG (low pressure side) 


SIGH 




Initial value of SIG (high pressure side) 


SL 




Slope in linear interpolation 


SLOPE 




Subroutine name 


SOL 




Stagnation entropy at left end (external) 


SOR 




Stagnation entropy at right port 


START 




Label for common block 


STEP 




Subroutine name 


T(I) 


T 


Temperature 


THETA(I) 




Angles in diffuser geometry 


TI 




Initial temperature in duct (low pres- 
sure side) 


TII 




Initial temperature in duct (high pres- 
sure side) 


TIM 


t 


Normalized time 


TMAX 




Cycle time 


TO 




Stagnation temperature 


TOOL 




Stagnation temperature in duct, left 
boundary 


TOCR 




Stagnation temperature in duct, right 
boundary 


TOL 




Stagnation temperature at left end 
(external ) 


TOP(I) 




Stagnation temperature at duct end 
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1 



:! 



"I 



I 

\ 






■f 




TOR 



U(I) 

X 

X(I) 

X1(I) 

X2(I) 

X3(I) 



Stagnation temperature at right end 
(external ) 

u Velocity 

Dummy variable 

X Spatial distance 

X interpolated along BEN 

X interpolated along BEM 

X interpolated along BES 
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PLCTTI.NG AND PRINTING PARAMETERS 
DATA TMAX/.3D0/ 

DATA AK,iM,L,NPR1NT/1.0,100, 1,02/ 
DATA RAN/1. 0,0. ,20. 0,0./ 

DATA NPLOT/100/ 
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